In [165]:
import numpy
import matplotlib
import numpy as np
import matplotlib.cm as cm
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
%matplotlib inline
from pylab import *
Liddriven Cavity:
Before exploring this example, note that we changed the basic code and removed loops. This really helps us to save the time. If you run this code with considering loop instead of indexes it will take approximately 15 hours to run. In this simulation we have a 2D fluid flow that is driven by a lid at the top which moves at a speed of. The cavity has four sides which three of them are solid and no-slip boundary condition is considered for walls and the top one which has moving velocity. If you remember in the file (illustration) we mention how to deal with moving or solid boundaries. This problem is two dimensional. We need to consider both advection and diffusion roles in Navier-stokes to ideally simulate this problem. In this example in contrast to previous examples the fluid has the ability to move.
Question:
How can we measure the strength of motion or strength of the fluid induced by moving the top lid? We have the answer thanks to the fluid mechanic course. The fluid could be Water or Gas with different kinematic velocities or densities. We at first need to non-dimensionalize the Navier-stokes equation. If we do the Reynolds number will appear which is defined as the inertia force over viscosity force. Higher Reynolds number implies higher inertia or lower viscosity resistance and vice versa.
What changes should be considered to shift our code from conduction to fully Navier-Stokes equation?
1- Define parameter “Reynolds” which you want to consider, it means that you definitely understand that you cannot directly use the parameters which you use in FDM or in original format of your equation. For example if you set the Reynolds number with adjusting the top velocity or kinematic viscosity in original equation here you are not allowed to consider any arbitrary value as same as you choose in FDM . The most important point is that you must make equal value of Reynolds number. You need to have Reynolds number for example 1000 then you chose the characteristic velocity as 0.1 (large value breaks low Mach number) and also you need to have relaxation time lower than 0.5 (Chapman-Enskog for Navier stokes suggests us). We extract the viscosity from Reynolds definition after adjusting characteristic velocity (lid velocity) and characteristic length (ny- the number of lattice nodes in y direction). By obtaining kinematic viscosity relaxation time could be obtained by satisfying the mentioned criteria.
2- We already set D2Q9 lattice velocities so we do not need to change streaming section however collision distribution function should be modified to account the advection roles. (Make attention how we change this part). 3- Comparisons were made by "Ghia.et al”. We used the u velocity along the horizontal line and V velocity along vertical line at the middle of the cavity .We got these data from “Gerris Flow Solver" site.
Your emission:
1- Try to change the Reynolds number to see what would be the vortex behavior. 2- Try to change the number of nodes in Lattice. What is the relation between lattice nodes and the length of the cavity? 3- Try to change the lid driven position (for example change the left side from no-slip to moving boundary). 4- Try to create cavity with unequal sides' length. Let’s the length of the cavity is twice of the width.
Equations: $$\frac{\partial \rho}{\partial t} + \overrightarrow{\nabla}\cdot(\rho\overrightarrow{u})=0 $$
$$\frac{\partial(\rho \overrightarrow{u})}{\partial t} + \overrightarrow{\nabla}\cdot[\rho\overline{\overline{u\otimes u}}] = -\overrightarrow{\nabla p} + \overrightarrow{\nabla}\cdot\overline{\overline{\tau}} $$and
Re= $\frac{u_0 \cdot ny}{\nu_{lbm}}$
Take care that you should necessarily make equal Reynolds number; It means that if you have a problem with Reynolds number 400 then it requires that:
$$Re_{lbm} = \frac{u_0 \cdot ny}{\nu_{lbm}} = \frac{u_{lid}\cdot H}{\nu} = 400 $$
In [166]:
nx=100 # the number of nodes in x direction lattice direction
ny=100 # the number of nodes in y direction lattice direction
Re=1000 # rayleigh number
uo=0.1
rhoo=1.0
visco=uo*ny/Re
D=3.0 # the dimension of the problem
omega=1./(0.5+(visco*D)) #Chapman-Enskog dt=1 and dx=1 relaxation flow
mstep=100000 # the number of time step
k=9 # k=0,1,2,3,4,5,6,8,9
In [167]:
w=numpy.ones(k) # witghting factor
cx=numpy.ones(k)
cy=numpy.ones(k)
rho=numpy.ones((ny+1,nx+1) ) # density matrix
u=numpy.ones((ny+1,nx+1) )
v=numpy.ones((ny+1,nx+1) )
f= numpy.ones((k, ny+1,nx+1)) # distribution function
feq= numpy.ones((k, ny+1,nx+1))
t1=numpy.ones((ny+1,nx+1) )
t2=numpy.ones((ny+1,nx+1) )
rhon=numpy.ones(nx+1)
In [168]:
y1= numpy.zeros(17)
u1= numpy.zeros(17)
y1[0]=-0.49933 ;u1[0]= -0.000882
y1[1]=-0.444335 ;u1[1]= -0.181701
y1[2]=-0.43629 ;u1[2]= -0.201989
y1[3]=-0.428914 ;u1[3]= -0.222276
y1[4]=-0.397406 ;u1[4]= -0.297251
y1[5]=-0.327052 ;u1[5]= -0.383699
y1[6]=-0.217948 ;u1[6]= -0.27788
y1[7]=-0.046595 ;u1[7]= -0.106804
y1[8]=0.001598 ;u1[8]= -0.060949
y1[9]=0.118733 ;u1[9]= 0.057217
y1[10]=0.235193 ;u1[10]= 0.186849
y1[11]=0.352315 ;u1[11]= 0.333239
y1[12]=0.45404 ;u1[12]= 0.466401
y1[13]=0.461386 ;u1[13]= 0.511382
y1[14]=0.469392 ;u1[14]= 0.574884
y1[15]=0.476719 ;u1[15]= 0.659554
y1[16]=0.5 ;u1[16]= 0.999118
In [169]:
x1= numpy.zeros(17)
v1= numpy.zeros(17)
x1[0]=-0.5 ;v1[0]=0.00069404
x1[1]=-0.43768 ;v1[1]=0.275621
x1[2]=-0.429602 ;v1[2]=0.290847
x1[3]=-0.421523 ;v1[3]=0.303994
x1[4]=-0.406521 ;v1[4]=0.326826
x1[5]=-0.343624 ;v1[5]=0.371038
x1[6]=-0.273803 ;v1[6]=0.330015
x1[7]=-0.265724 ;v1[7]=0.32307
x1[8]=-0.000289 ;v1[8]=0.0252893
x1[9]=0.304962 ;v1[9]=-0.318994
x1[10]=0.359781 ;v1[10]=-0.427191
x1[11]=0.40652 ;v1[11]=-0.515279
x1[12]=0.445182 ;v1[12]=-0.392034
x1[13]=0.45326 ;v1[13]=-0.336623
x1[14]=0.461339 ;v1[14]=-0.277749
x1[15]=0.46884 ;v1[15]=-0.214023
x1[16]=0.5 ;v1[16]=-6.20706e-17
In [170]:
strf=numpy.ones((ny+1,nx+1) ) # streamfunction
##================ Initial boundary condition
w[0]=4./9 #4./9
w[1:5]=1./9.
w[5:9]=1./36.
cx[0]=0.0
cx[1]=1.0
cx[2]=0.0
cx[3]=-1.0
cx[4]=0.0
cx[5]=1.0
cx[6]=-1.0
cx[7]=-1.0
cx[8]=1.0
cy[0]=0.0
cy[1]=0.0
cy[2]=1.0
cy[3]=0.0
cy[4]=-1.0
cy[5]=1.0
cy[6]=1.0
cy[7]=-1.0
cy[8]=-1.0
In [171]:
##================== Initial value
rho[0:ny+1,0:nx+1]=rhoo
v[0:ny+1,0:nx+1]=0.0
u[0:ny+1,0:nx+1]=0.0
u[ny,0:nx+1]=uo
for i in range(nx+1):
for j in range(ny+1):
for l in range (k): #k=0,1,2,3,4
f[l,j,i]=w[l]*rho[j,i]
In [172]:
## Main loop : comprised two parts :collision and streaming
##=====================
for n in range(mstep) :
#print n ,u[nx/2,ny/2]
if (n%10000==0):
print n
# collision process fluid flow
t1[0:ny+1,0:nx+1]=u[0:ny+1,0:nx+1]*u[0:ny+1,0:nx+1]+v[0:ny+1,0:nx+1]*v[0:ny+1,0:nx+1]
for l in range (k):
t2[0:ny+1,0:nx+1]=cx[l]*u[0:ny+1,0:nx+1]+cy[l]*v[0:ny+1,0:nx+1]
feq[l,0:ny+1,0:nx+1]=w[l]*rho[0:ny+1,0:nx+1]*( 1.+ 3.*t2[0:ny+1,0:nx+1]+4.5*t2[0:ny+1,0:nx+1]**2-1.5*t1[0:ny+1,0:nx+1] )
f[l,0:ny+1,0:nx+1]=(1.-omega)*f[l,0:ny+1,0:nx+1]+omega*feq[l,0:ny+1,0:nx+1]
f[1,0:ny+1,nx:0:-1]=f[1,0:ny+1,nx-1::-1]
f[3,0:ny+1,0:nx]=f[3,0:ny+1,1:nx+1]
f[2,ny:0:-1,0:nx+1]=f[2,ny-1::-1,0:nx+1]
f[5,ny:0:-1,nx:0:-1]=f[5,ny-1::-1,nx-1::-1]
f[6,ny:0:-1,0:nx]=f[6,ny-1::-1,1:nx+1]
f[4,0:ny,0:nx+1]=f[4,1:ny+1,0:nx+1]
f[8,0:ny,nx:0:-1]=f[8,1:ny+1,nx-1::-1]
f[7,0:ny,0:nx]=f[7,1:ny+1,1:nx+1]
## =============================
#left Boundary- Bounce back
f[1,0:ny+1,0]=f[3,0:ny+1,0]
f[5,0:ny+1,0]=f[7,0:ny+1,0]
f[8,0:ny+1,0]=f[6,0:ny+1,0]
#right Boundary-bounce back
f[3,0:ny+1,nx]=f[1,0:ny+1,nx]
f[7,0:ny+1,nx]=f[5,0:ny+1,nx]
f[6,0:ny+1,nx]=f[8,0:ny+1,nx]
# top Boundary -bounce back
rhon[0:nx+1]=f[0,ny,0:nx+1]+f[1,ny,0:nx+1]+f[3,ny,0:nx+1]+2*(f[2,ny,0:nx+1]+f[6,ny,0:nx+1]+f[5,ny,0:nx+1])
f[4,ny,0:nx+1]=f[2,ny,0:nx+1]
f[7,ny,0:nx+1]=f[5,ny,0:nx+1]-(rhon[0:nx+1]*uo/6.0)
f[8,ny,0:nx+1]=f[6,ny,0:nx+1]+(rhon[0:nx+1]*uo/6.0)
# bottom Boundary -bunce back
f[2,0,0:nx+1]=f[4,0,0:nx+1]
f[5,0,0:nx+1]=f[7,0,0:nx+1]
f[6,0,0:nx+1]=f[8,0,0:nx+1]
for i in range(nx+1):
for j in range(ny+1):
sumr=0.0
for l in range (k):
sumr=sumr+f[l,j,i]
rho[j,i]=sumr
for i in range(1,nx):
for j in range(1,ny):
usum=0.0
vsum=0.0
for l in range (k):
usum=usum+f[l,j,i]*cx[l]
vsum=vsum+f[l,j,i]*cy[l]
u[j,i]=usum/rho[j,i]
v[j,i]=vsum/rho[j,i]
u[:,0]=0.
u[:,nx]=0.
u[0,:]=0.
u[ny,:]=uo
v[:,0]=0.
v[:,nx]=0.
v[0,:]=0.
v[ny,:]=0.
In [173]:
strf[0,0]=0.0
for i in range(nx+1):
rhoav=0.5*(rho[0,i-1]+rho[0,i])
if(i>0) :
strf[0,i]=strf[0,i-1]-rhoav*0.5*(v[0,i-1]+v[0,i] )
for j in range(1,ny+1):
rhom=0.5*(rho[j,i]+rho[j-1,i])
strf[j,i]=strf[j-1,i]+rhom*0.5*(u[j,i]+u[j-1,i])
x=numpy.linspace(0,nx,nx+1)
y=numpy.linspace(0,ny,ny+1)
In [174]:
mx, my = numpy.meshgrid(x,y)
plt.figure(figsize=(8,5))
plt.contourf(x,y,strf,20)
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.colorbar();
In [190]:
mx, my = numpy.meshgrid(x,y)
plt.figure(figsize=(8,5))
plt.contour(x,y,strf,900)
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.colorbar();
In [175]:
mx, my = numpy.meshgrid(x,y)
plt.figure(figsize=(8,5))
plt.contourf(x,y,u,20)
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.colorbar();
In [176]:
mx, my = numpy.meshgrid(x,y)
plt.figure(figsize=(8,5))
plt.contourf(x,y,v,20)
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.colorbar();
In [177]:
x1=(x1*ny)+(ny/2)
y1=(y1*ny)+(ny/2)
vv=v[ny/2,0:nx+1]/uo
uu=u[0:ny+1,nx/2]/uo
plt.plot(x1,v1, 'r*',label=' Ghia et.al');
plt.plot(x,vv, 'b-',label=' LBM ');
plt.legend();
In [178]:
plt.plot(y1,u1, 'r*',label=' Ghia et.al');
plt.plot(y,uu, 'b-',label=' LBM ');
plt.legend();